Theory of Structural Response to Macroscopic Electric Fields in Ferroelectric Systems 
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We have developed and implemented a formalism for computing the structural response of a 
periodic insulating system to a homogeneous static electric field within density-functional perturba- 
tion theory (DFPT). We consider the thermodynamic potentials E(R,,n,£) and i ? (R, ?7,P), whose 
minimization with respect to the internal structural parameters R and unit cell strain rj yields the 
equilibrium structure at fixed electric field £ and polarization P, respectively. First-order expansion 
of -E(R, ri,£) in £ leads to a useful approximation in which R(P) and j?(P) can be obtained by 
simply minimizing the zero-field internal energy with respect to structural coordinates subject to 
the constraint of a fixed spontaneous polarization P. To facilitate this minimization, we formulate 
a modified DFPT scheme such that the computed derivatives of the polarization are consistent with 
the discretized form of the Berry-phase expression. We then describe the application of this approach 
to several problems associated with bulk and short-period superlattice structures of ferroelectric ma- 
terials such as BaTi03 and PbTi03. These include the effects of compositionally broken inversion 
symmetry, the equilibrium structure for high values of polarization, field-induced structural phase 
transitions, and the lattice contributions to the linear and the non-linear dielectric constants. 



PACS numbers: PACS numbers: 77.22. Ch, 77.65. Bn, 77.84.Dy, 71.15.-m 
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I. INTRODUCTION 

As the usefulness of density-functional theory (DFT) 
for the study of dielectric materials is now well estab- 
lished, one might imagine that calculations of crystalline 
insulators in the presence of a homogeneous macroscopic 
electric field should be routine. On the contrary, the 
presence o£|-an electric field introduces several severe 
difficulties .Era The electric potential acquires a term that 
is linear in the spatial coordinates, thus violating the pe- 
riodicity condition underlying Bloch's theorem and act- 
ing as a singular perturbation on the electronic eigen- 
states. Moreover, in principle there is no longer a well- 
defined ground state for the electrons in a solid in a 
macroscopic electric field because the energy of the sys- 
tem can be always lowered by transferring electrons from 
the valence band in one spatial region to the conduction 
band in a distant region. 

One way around these difficulties is to make 
of density-functional perturbation theory (DFPT) 
which provides a framework for calculating the pertur- 
bative response to infinitesimal electric fields (as well as 
to atomic displacements and strains). DFPT has been 
widely adopted for many studies of the dielectric and 
piezoelectric properties and dynamic effective charges of 
dielectric materials. However, being a perturbative ap- 
proach, the method is not capable of treating a finite 
electric field directly. 

A more direct attack on the-finite-field problem was 
made by Nunes and Vanderbilt, El who showed that a real- 
space Wannier-function representation could be used to 
represent the electronic system in the presence of a fi- 
nite electric field.EI In this approach, one minimizes a 
total-energy functional of a set of field-dependent Wan- 
nier functions for a periodic system at fixed electric field. 
Alternatively, the minimization can also be performed at 



fixed polarization via a Legendre transformation of the 
energy functional with the electric field treated as a La- 
grange multiplier. The approach was implemented in the 
DFT context by Fernandez et al.p but proved too cum- 
bersome to find widespread utility. 

In the present paper we propose a new scheme for the 
treatment of a dielectric system in a static homogeneous 
electric field. Our scheme is based on a low-order trunca- 
tion of the DFPT perturbative expansion in electric field, 
and the use of this truncated expansion to extrapolate to 
finite electric field. A key feature of our approach is that, 
while we keep only low orders in the expansion in elec- 
tric field, we effectively keep all orders of expansion in 
the structural degrees of freedom. We demonstrate that 
even a first-order truncation of the electric-field pertur- 
bation provides a very useful and practical scheme. In 
this context the electric field simply couples to the zero- 
field polarization, so that the latter plays a central role 
in our formulation. In fact, it is rather natural to formu- 
late our approach in terms of a constrained minimization 
procedure in which the DFT energy functional is mini- 
mized over all structural degrees of freedom subject to 
a constraint on the value of the polarization. This al- 
lows a two-step approach in which one first maps out 
the energy surface as a function of polarization in the 
DFT framework, and then uses this energy surface, aug- 
mented by the coupling to the electric field, to obtain the 
ground-state structure in the presence of the field. We 
will show that essentially no additional approximations 
are needed beyond the first-order truncation of the free- 
energy expansion in electric field. We will also show that 
the methodology can be extended to second order (or, in 
fact, to any desired order) in the electric field. 

Before proceeding, we should acknowledge an addi- 
tional theoretical subtlety associated with the correct 
choice of exchange-correlation functional in the electric- 
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field problem. Gonze, Ghosez and GodbyB (GGG) ar- 
gued that the exact Kohn-Sham exchange-correlation 
energy functional should have a dependence on the 
macroscopic polarization and formulated a "density- 
polarization functional theory" in which there is gener- 
ally an exchange-cor.celation contribution to the Kohn- 
Sham electric fieldE While this was an important 
formal d^elopment that subsequently received much 
attentionBBIilJ it has not yet led to an improved prac- 
tical exchange-correlation functional. We thus restrict 
ourselves here to the usual LDA exchange-correlation 
functional where, because of the locality of the central 
approximation, the subtleties identified by GGG do not 
arise. 

This paper is organized as follows. In the next section, 
we present our formalism for computing the structural re- 
sponse of an insulating system to an electric field. Some 



details of the implementation are presented in Sec. Ill 



including details of our minimization procedure, a discus- 
sion of modifications that we made to the DFPT proce- 
dure to achieve compatibility with the discretized Berry- 
phase polarization formula, and a description of the tech- 
nical details of the ab-initio pseudopotential calculations. 
Then, in Sec. IV, we pre sent s everal sample applications 
of our method. In Sec. [VA, we show that it provides 



an alternative approach to the study of short-period fer- 
roelectric piperlatti ce st ructures with broken inversion 
symmetryO In Sec. IV B we present a study of the depen- 
dence of the internal struc tural parameters of BaTi03 on 
polarization. In Sec. IV C| we show that our method pro- 
vides a straightforward way of computing the dielectric 
susceptibilities and piezoelectric coefficients as functions 
of the electric field, thus allowing an estimation of the 
non-linear dielectric and piezoel ectric response in a ferro- 
electric system. Finally, in Sec. IV D| , we consider a case 
in which a full three-dimensional treatment of the polar- 
ization and the structural distortions is needed. Specifi- 
cally, we model the polarization rotation and structural 
phase transitions induced by the application of a macro- 
scopic electric field to a model ferroelectric system_and 
relate our results to recent experiments in PZN-PT.Ej Fi- 
nally in Sec. [V], we summarize our work and discuss the 
prospects for future applications of our approach. 



II. FORMALISM 

Our goal is to investigate the effect of a homogeneous 
static electric field on the structure and polarization of 
polar insulators, including systems with a nonzero spon- 
taneous polarization. In addition to an efficient approach 
for computation, we also aim towards a formulation that 
readily allows an intuitive understanding of the effects 
of the field. As will become clear below, this will lead 
us to a formulation in which the polarization plays an 
especially prominent role. 



A. Case of a single minimum 

Let £ be the macroscopic electric field, R the atomic 
coordinates, and r\ the lattice strain, and assume that 
the total energy per unit cell E(R, rj,£) has a single lo- 
cal minimum of interest in the (R, rf) space for given £ . 
(This restriction is normally appropriate for a paraelec- 
tric material, but not for a ferroelectric one. The ex- 
istence of multiple local minima in the latter case calls 
for a more careful discussion, which is deferred to the 
following subsection.) We let 



E(£) = min E(R, 77, £) 

H.,7] 



(1) 



and denote the location of the minimum by R cq (£) 
and 7y cq (£). The polarization P(R, ?7,£), the thermody- 
namic conjugate of £, can then be obtained from the 
expression P(R, rj,£) = — [dE(R, rj,£)/ d£]n. v and P(£) 
obtained by evaluating P(R oq (£), 7y eq (£), £), or equiva- 
lently, -dE{£)/d£. 

We can recast this minimization into a form in which 
the polarization is more central. Viewing i?(R, rj,£) as 
a thermodynamic potential that minimizes to equilib- 
rium values of R and 77 at fixed £ leads naturally via 
a Legendre transformation to a thermodynamic poten- 
tial -F(R, ?/, P) that minimizes to equilibrium values of 
R and 77 at fixed P: 

F(R,77,P) = min[-E'(R, 17, A) + A • P] 

= £(R,77,A(R,77,P)) + A(R, ?i ,P).P(2) 

with A(R, 77, P) being the value at the minimum. This 
is equivalent to P(R, 77, A) = P, that is, A(R, 77, P) is 
the value of the macroscopic field necessary to produce 
polarization P at given R and 77. 
We then define the function 



F(P) = minF(R,?y,P) 

R,?7 

= ^(R cq (P),77 cq (P),P) 



(3) 



with R oq (P) and ?7 cq (P) being the values at the mini- 
mum. These structural parameters R oq (P) and 77 oq (P) 
are in fact equal to R eq (£) and ?7 eq (£), the structural 
parameters defined by minimizing E(R, 77, £) at the cor- 
responding fixed £ = A(R cq (P), 77 cq (P), P). The polar- 
ization at this extremum is, as expected, 

P(R(£), ri(£), £) = 
P(R cq (P),77 eq (P),A(R cq (P),77 cq (P),P)) = P. (4) 

Finally, we re-express the original minimization as 

E(£) = min[F(P) - £ • P] . (5) 

In this expression, the electric field £ appears only in the 
term — £ ■ P, and thus the effects of the field can be com- 
pletely understood by investigating the ^-independent 
free energy landscape F(P). 
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B. Case of multiple stationary points 

In many cases of interest, the function E(R.,r],£) has 
several local minima, and the essential physics of the 
problem is to map out the competition between these 
minima. For example, in a tetragonal ferroelectric like 
PbTiC>3, there are six degenerate minima of this func- 
tion at £ = 0, and the application of a nonzero £ breaks 
the symmetry and establishes one dominant domain ori- 
entation of the polarization. However, it may also be of 
interest to follow the behavior of the other local minima, 
corresponding to metastable states, as well as other sta- 
tionary points of this energy surface. For example, saddle 
points and local maxima of E(R, 77, £ ) can correspond to 
stable states for fixed P. 

In such cases, it is straightforward to generalize the 
previous discussion by associating a label (n) with each 
stationary point of interest. Thus the location of the sta- 
tionary point is denoted by PJ™)(£) and and 
= £(P>)(£),77( n ) (£),£) is the energy at the 
stationary point. The arguments of the previous sub- 
section carry over much as before. The discussion fol- 
lowing Eq. ([5]) is modified by noting that the minimiza- 
tion of F(R, 77, P) with respect to R and 77 at fixed P in 
Eq. (||) will always be associated with one of the station- 
ary points of E(R, 77, £ ) with respect to R and 77 at the 
corresponding fixed £; that is, R cq (P) = R'™)(£) and 
7? cq (P) = rf n >{£) for some n. 

Finally, defining the global minimum E{£) = 
mm n E^ n \£), it is easy to see that Eq. (0) holds as be- 
fore. 



C. Truncation of the expansion 

The central quantities appearing in the preceding sub- 
sections are the energy E(R, rj, £) and the polarization 
P(R, 77, £) in a given electric field £. Unfortunately, there 
is as yet no rigorous formulation of DFT for the case 
of finite non-zero £ . However, electric-field derivatives 
of arbitrary order may be computed by the methods of 
density-functional perturbation theory. We thus expand 
in £ around £ = 0: 



E(R )V ,£) = E(R,r,,£ = 0) + J2^ dE{ f^' £) 



d 2 E(R,7 l7 £) 



a/3 



d£ a d£ 







(6) 



£=0 



Carried to all orders in £ , this expansion is exact. How- 
ever, for sufficiently small fields we can make the approx- 
imation of truncating this sum to define Bi(R, 77, £) as 
the sum of the first i+1 terms in Eq. (^|), and Pj(R, rj, £ ) 
as — (<£Ej(R, rj, £)/d£)\ R . Note that this truncation is 
only in powers of £, and that the dependence on R and 
77 is preserved to all orders. 



For many systems, it is already of interest to consider 
the simplest case i = 1, where 

£i(R, 77, £) = E(R, 77, 0) - £ ■ P(R, 77, 0) , (7) 



P!(R,77,£)=P(R,77,0). 
At this order, the resulting expression 



(8) 



F(P) = min [E(R, r), 0) + A • (P - P(R, 77, 0))] (9) 

R,?7,A 

can also be interpreted as one in which A appears sim- 
ply as a Lagrange multiplier implementing the constraint 
P(R, ?7, 0) = P in the set of equations that minimize 
E(R, rj, 0) over R and 77, i.e., 



<9£:(R,77,0) ^ dPp(TL,i},0) 



dR,., 



^ dRi, 

(3 



aE(R,77,0) ^^(R,7?,0) 



P(R,77,0)=P 



P 



A/3 = 



Xb = 



(10) 



Details of the calculation are described in the next sec- 
tion, and all of the results reported in the following sec- 
tions are obtained using the i — 1 expressions. 

Generalization of the formalism to order i = 2 is pro- 
vided in Appendix 



D. Relation to method of Fu and Cohen 

In the remaining part of this section we discuss-an ear- 
lier approach proposed by Fu and Cohen (FC)Jl3 These 
authors carried out a first-principles investigation of the 
mechanism of rotation of the polarization in BaTi03 by 
an applied electric field. Their approach is similar to 
our i — 1 case, but involves an additional approximation 
which we will describe by expressing their procedure in 
the notation established in this section. 

The first step in their approach is the same as our 
7 = 1 case: to approximate E(TL, 77, £) by £i(R, r],£). 
The next step is to perform a constrained minimization, 
computing 



U(Q) 



min £'(R, 77,0) 

Q(R, 7? )=Q 



(11) 



where the constraint is not on the polarization, but on 
Q(R, 77) where Q(R, 77) is the Ti displacement relative 
to the average position of the other atoms in the unit 
cell. R(Q) and 7/(Q) will be defined as the values of the 
atomic coordinates and strain at the minimum. 

The equilibrium energy E-pc(£), structure 
R(Qmin), »7(Qmin) and polarization Pfc(£) = 
P(R(Q min ), 77(Q min ), 0) are then obtained by the 
minimization 



£ F c(£)=min{f/(Q) 
Q 



£.P(R(Q),7j(Q),0)}. (12) 
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However, the results will in general not be equal to those 
obtained with our i = 1 expression 

E{£) = min{ min ECR, 77, 0) - £ ■ P} . (13) 

P P(R,?j,0)=P 

The reason is that 
E{£) = min{£;(R,77,0)-f -P(R,77,0)} 

R,77 

= min{ min {ECR, 77, 0) - £ ■ P(R, n, 0)}} 

Q Q(R,77)=Q 

< min{ {ECR, 77, 0) - £ ■ P(R, 77, 0)}| R(Q)i>j(Q) } 



^Fc(f)- 



(14) 



The point is that the coordinates R(Q) and r](Q) that 
minimize E alone at fixed Q are generally not the same 
as those that would minimize the combination (E — S- P) 
at fixed Q. At nonzero £ , equality will only be obtained 
under very special circumstances, for example, if the sur- 
faces of constant P in R, 77 space coincide with the sur- 
faces of constant Q, at least for the relatively low energy 
structures. More specifically, the polarization should be 
a function only of the Ti displacement relative to all other 
atoms independent of the detailed arrangement of those 
atoms and of the strain. This is not unreasonable for 
small enough fields in BaTiC>3, in which the soft mode 
is almost a pure Ti displacement and is well isolated in 
energy from other polar modes. However, as we will see 
in the following discussion, this proportionality is never 
quite exact even for small fields, and the discrepancies 
grow rapidly as the fields get larger. 



III. METHODOLOGY 
A. Minimization procedure 

We now describe in detail how the minimization of 
E(R, 77,0) respecting the constraint P(R, 77, 0) = P is 
implemented. Eq. (|9|) shows that this constraint can 
be imposed by a Lagrange multiplier A. Therefore, we 
need to solve the stationary-value problem described by 
Eqs. ©. 

Suppose we make a trial guess of the initial coordinates 
Ro and strains 770 for the desired structure. The energy 
E(R, rj, 0) can be expanded up to second order in 5R = 
R — R and 8rj = 77 — 770 as 

£(R,77,0) = E(R o , Vo ,0) +J2(-F ta )SR ia 

ict 

+ 9 Y c H» S Vn$Vi> + 'Yia SR *'* ST ll* ( 15 ) 



where Fi a are the Hcllmann-Feynman forces, cr M are the 
stresses in Voigt notation, K % K are the force-constant 



matrix elements, c M „ are the elastic constants, and 7; QM 
are the coupling parameters between the internal coor- 
dinates and strains. The corresponding variation in the 
polarization P(R, 77, 0) is 

P Q (R,77,0) = P Q (R ,77o,0) + ^z; / 3«5i 1 , J/ 3 + ^e Q 



where Z 1 * = dP a /dRip and e afl — dP a /drj^ are re- 
spectively the dynamic effective charge and piezoelectric 
tensors. 

Eqs. ( |Io|) lead to the linear system of equations 



7 c e 
Z* e 



SR\ 




f F 


Sri 


-1 


a 


A / 







(17) 



for 5R, Sjj and A, where V on the right-hand side denotes 
the difference between the initial and target values of P. 
At each step of the minimization, we compute SR and 6r}, 
and obtain the new coordinates and strains via RJj ew = 
Ro + SR and t/q™ = 770 + ^77. Then R ncw and 77 ncw are 
chosen as the new "trial" coordinates and strains. This 
is repeated until convergence is achieved. 

For a practical implementation of this procedure we 
use density-functional perturbation theory, which allows 
us to compute the coefficients K l £a and Z l a p efficiently. 
The forces F and the stresses a are calculated bsz-.the 
Hellmann-Feynman theorem with Pulay correctionalj for 
the stresses. However, the computation of the remaining 
quantities in Eq. (17), involving derivatives with respect 
to strain, requires two additional comments. First, the 
DFPT calculation of 7, c and e is not yet implemen ted in 
the current version of the ABINIT package (see Sec. [II C| ) 
and the finite-difference calculation of these quantities 
would be exceedingly tedious. However, we will show 
that an alternative indirect minimization can be carried 
out by means of a Devonshire -type Ha milton ian.Ej De- 
tails will be given in Section [V C and [VP . Second, 
the most efficient way to compute V is with a discretized 
Berry-phase expression. However, the dependence on R 
and 77 of the resulting polarization corresponds to the 
DFPT derivatives exactly only in the limit of a dense k- 
point sampling mesh. This issue is discussed and resolved 
in detail in the next subsection. 

Before concluding this subsection, we note that the 
higher-order formalism can be implemented in an analo- 
gous way. However, additional energy derivatives would 
be needed. The details of the treatment for i = 2 are 
presented in Appendix For the following study, we 
will restrict ourselves to the first-order case described by 
Eq. 



B. DFPT computation of derivatives of the 
discretized Berry-phase polarization 

In the implementation of the minimization procedure 
(Eq. |j"7| ), a practical problem arises in connection with 
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the calculation of the dynamical effective charges Z* and 
polarization P. By definition, they should be related by 



J a/3 



V 



dP a 
dR 



iff 



(18) 



where a and j3 are Cartesian directions, i is the index for 
the atom, and V represents the unit cell volume. 

However, when the discretized Berry-phase expression 
is used to compute P and the DFPT expression is used 
to compute Z* on the same k-point mesh, Eq. ( fL8| ) is not 
satisfied exactly. The discrepancy vanishes in the limit 
of a dense k-point mesh, but in a practical calculation, 
which must use a finite mesh, it will result in an incon- 
sistency in the equations for.-the minimization. 

In the Berry-phase theorytUj the polarization is 



dv xc (r)/dRi a are, respectively, the first-order derivatives 
of the external potential and the exchange-corcelation 
potential with respect to a q = displacement Jla The 
third expression (omitted here) includes both types of 
first-order wavefunction changes, and is stationary with 
respect to small errors in the first-order wavefunctions. 

In DFPT, first-order changes IV^ 1 ^) in the wavefunc- 
tions with respect to a perturbation can be computed 
as self-consistent solutions of the first-order Sternheimer 
equations^ 



P C {H 



)Pc\4>. 



(i) 



-Pc# (1) |V> 



(0)\ 



(23) 



subject to a "parallel transport" gauge constraints 

(^°W>=0, (24) 



ife 
(2tt) 3 



BZ 



occ 

E 



«mk )dk 



(19) 



where f — 2 for spin degeneracy. P BP is computed in 
practice using a discretized formula which, for the case 
of isolated bands, takes the form 



Pn 



/ dk±y^ Iraki JJ (w mk |um,k+b 



(2 



ke5(kj_) 



(20) 



where the integration over the 2D plane perpendicular 
to direction a is replaced in practice by a summation over 
a 2D mesh. The product runs over a string <S(k_i_) of k- 
points running parallel to direction a at a given kj_. b 
is the separation between neighboring points along the 
string and / = 2 is the spin degeneracy factor. The 
composite-band formulation is presented in Appendix 

In DFPTE3 there arc three equivalent expressions for 
the Z* tensor. The first is the change in the polariza- 
tion due to the first-order change in the wavefunctions 
resulting from an atomic displacement: 



Z 



a/3 



ifeV 
(2tt) 3 



BZ 



occ 

E 



dR,,, 



du ri 



dk 



dk, 



(21) 



where du m k/dRi a is the first-order change of the wave- 
function due to the perturbation by displacing an atom 
belonging to the ith sublattice along the a axis. Alterna- 
tively, Z* can be computed as the derivative of the force 
along direction a on an atom in the ith sublattice with 
respect to an electric field along direction /?, 



Zap — 2 



V 



OCC 



2tT Jbz 
dv xc (r) dn(r) 
dRi a d£f3 



f(u. 



mk 



dRia 



d£ 6 



■dr 



rfk 



(22) 



where du^/ d£p is the first-order change of the wave- 
functions due to the electric field, and dv CK t/dRi a and 



where ifW is the first-order change in H and P c is the 
projection operator onto the subspace of the conduction 
bands, and n and m run only over the valence bands. 

In the case of the electric-field perturbation with field 
in Cartesian direction a, the Sternheimer equation takes 
the forirMN 



P C (H - e mk )P c 



du ri 



= -P c H^\u mk ) 



where 



ff« = -i 



dk a 



dvn 
d£ a 



dv xc (r) 
d£ a 



(25) 



(26) 



As input to this equation, we need the quantity 
du m k/dk a appearing on the right-hand side. This is ob- 
tained by solving a second Sternheimer equation 



Pc(H 



^77? k 



du Tl 



dk a 



Pr 



dk a 



(27) 



where = |(— iV + k) 2 + vks an d thus dH^/dka = 
— iVq, + k a . The presence of the operator id/dk a in 
Eq. (|2^) is a unique feature that arises from the cou- 
pling £ ■ P between the macroscopic electric field and the 
polarization. 

Now we come to the main point of this subsection, 
which is that the derivatives of P calculated from DFPT 
are not exactly the derivatives of the discretized Berry- 
phase expression for P in practical calculations. In par- 
ticular, for a given k-point sampling, the Z* computed by 
solving for d\u m ^)Jdk a from Eq. (|27| ), and then comput- 
ing Z* via Eq. ( Pl| ) or via Eqs. (|25|) and (|2^), is not ex- 
actly equal to the Z* computed from finite differences of 
the Berry-phase polarization in Eq. (|20|). The numerical 
discrepancy between the Berry-phase polarization and 
the Z* computed from the DFPT affects the applica- 
tion of t he co nstrained minimization scheme proposed in 
Section III A because it introduces an inconsistency into 
the linear system in Eq. (p"7j). 

Our cure for this problem is to modify the algorithm 
by which d\u m -^)/dk a is calculated within the DFPT 
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framework. Instead of solving the Sternheimer equa- 
tion (p7|), we calculate d\u m k) / dk a from finite differences 
of the ground state Bloch wavefunctions |u m k) between 
the neighboring k-points along the a-direction. This ap- 
proach corresponds to the "perturbation expansion after 
discretization" formalism discussed by Nunes and Gonzc 
in Ref. | 

We illustrate our approach here for the case of a sin- 
gle band in one dimension. The generalization to the 
three-dimensional multi-band case is postponed to Ap- 
pendix |b| 

For a single band in one dimension, the Berry-phase 
polarization is 



P = — — 'S~^ Im ln< 

9tt 



2tt 



Uk\Uk+b) 



(28) 



and its first-order variation reflecting a first-order change 
in the wavefunction \5v,k) is 

~(5u k \u k+b ) 



SP 



2tt 



5> 



(u k \Su k -i 



9ir 



2tt 
feb 



{Uk\Uk+b) 

(6u k \u k+b ) 



Uk\Uk+b) 

y^Ke{8u k \v k ) 



[Uk\Uk+b) 

(6u k \u k _ b ) 
(uk\u k -b) 



where 



lVk) = Yb 



\Uk+b) 



K-b) 



(u k \u k+b ) (u k \u k - b ) 



(29) 



(30) 



is understood to be a finite-difference approximation to 
id\u k ) / dk. Note that Eqs. ( p9[pfj| ) are manifestly gauge- 
independent in the sense of being independent of the 
choice of phases for the \u k ). 

In the three-dimensional multi-band case, we just need 
to replace v k of Eq. ( |30|) by its generalization v m k,a rep- 
resenting durnk/dka as discussed in Appendix B; w r „k,a 
is gauge-independent in the more general sense of be- 
ing invariant with respect to a unitary rotation among 
occupied bands on neighboring k-points. This u m k,a can 
then be substituted for du m ^/dk a in Eq. ( ]2l| ) to compute 
Z* . Or equivalently, it can be inserted into Eq. (£5|) to 
compute du m k/d£ a , which in turn can be substituted 
into Eq. ( p2f ) to compute Z*. In either case, we are 
guaranteed to obtain the same values of Z* as would 
be derived from a series of finite-difference calculations 
of polarization vs. atomic displacement using the same 
k-point set. This is because Eqs. (|2^-p0|) are derived 
directly from the Berry-phase polarization expression of 
Eq. ( |l9| ) using the same k mesh. Moreover, because the 
Berry-phase polarization (including ionic contributions) 
is independent of origin, it also follows that the acoustic 
sum ruleEj Yli = <W on the components of the dy- 
namic effective charges will be satisfied exactly, which is 
not the case in conventional linear-response calculations 
otZ*. 



C. Computational Details 

We carried out— all the ab initio calculations using the 
ABINIT package,Ej in which we have implemented the 
above algorithm. ABINIT uses a plane-wave basis and 
provides multiple norm-conserving (NC) and extended 
NC pseudopotentials. The discretized formula for the 
wavefunction derivatives with respect to the wavevectors, 



Eq. (B14), is introduced in a new subroutine dudk.f, a 



key ingredient that allows us to carry out the constrained- 
polarizati on mi nimization scheme. 

In Sec. IV A, in order to, construct the pseudopoten- 
tials for the virtual atomstil that enter the heterovalent 
system Ba(Ti-<5,Ti,Ti+6)0 3 , we utilize the FHI atomic 
codeES that generates Troullier-Martin separable norm- 
conserving pseudopotentials. However, the FHI pseu- 
dopotential generation scheme only allows one projector 
within each angular momentum channel, thus preventing 
the inclusion of the 3s and 3p states, in addition to 3d 
and 4s states, in the valence for the Ti pseudopotential. 
(The same problem occurs for the 5s and 5p states for 
the Ba atom.) We generate the pseudopotential in ion- 
ized configurations 3s 2 3p 6 3<i 2 4s for Ti and 5s 2 5p 6 6s° for 
Ba. We used the exchange-correlation energy functional 
in the Ceperley- Aldeio form with Perdew and WangE£l 
parameterization. 

The studies described in Sections [VB to IV D 



have been performed with the highly transferable ex- 
tended-.norm-conserving pseudopotentials proposed by 
Teter.Ej A Perdew-ZungerEZI parameterized Ceperley- 
Alder exchange-correlation functional was used. These 
pseudopotentials include the Pb 5d, 6s and 6p, the Ba 
5s, 5p and 6s, the Ti 3s, 3p, Ad and 4s, and the O 2s and 
2p electrons in the valence states. 

We have used an energy cutoff of 35 Ha throughout. 
The integrals over the Brillouin zone have been replaced 
by a sum over a 4 x 4 x 4 k-point mesh. Both the k- 
point sampling and the energy cutoff have been tested 
for good convergence of the phonon eigenvalue and eigen- 
vector properties. We use the same k-point mesh for 
the Berry-phase calculations. Convergence of the relax- 
ations requires the Hellmann-Feynman forces to be less 
than 0.02 eV/A. (In the constrained minimization pro- 
cedure described by Eq. (|l7|), the forces that are tested 
for convergence are the ones after projection onto the 
constant-P subspace.) 



IV. SAMPLE APPLICATIONS 

In this section, we illustrate the theory within the first- 
order (i = 1) formalism (see Section II) by applying it to 
a series of problems involving ferroelectric, dielectric and 
piezoelectric properties. In particular, we emphasize that 
the main purpose of these calculations is to exhibit and 
understand the nonlinearity in the structural response of 
the ferroelectric systems to an electric field. Such studies 
have not previously been widely pursued. We have used 



7 




A. Inversion symmetry-breaking system 

In a conventional ABO3 perovskite such as BaTi03, 
the cubic symmetry of the high-temperature paraelectric 
phase is spontaneously broken at the transition to the 
ferroelectric phase. The atomic displacements that oc- 
cur in the ferroelectric phase give rise to an associated 
lattice strain, and the ferroelectric state is characterized 
by a switchablc polarization because of the occurrence 
of degenerate energy minima that are connected by the 
broken symmetry operations. 

Recently, using DFT-|total-energy methods, we (Sai, 
Meyer and Vanderbilttil) studied a new class of cubic 
perovskite compounds in which the composition is mod- 
ulated in a cyclic sequence of three layers on the A 
site (i.e., (AA 1 ' A")BO^ structures) or on the B site (i.e., 
A(BB' B")03 structures). The inversion symmetry that 
was present in the high-symmetry cubic structure is now 
permanently broken in these materials by the alternating 
compositions ordered along the lattice growth direction. 
This gives rise to important qualitative differences in the 
energetic behavior of these compounds relative to sim- 
ple ABO3 perovskites. Most interestingly, it was shown 
that by using heterovalent compositional substitutions, 
the strength of the breaking of the inversion symmetry 
could be tuned through an enormous range, suggesting 
that such systems could be very promising candidates for 
new materials with largB-piezoelectric and other dielec- 
tric response properties. til 

Such compositionally modulated structures were stud- 
ied within a model system Ba(Ti— <5,Ti,Ti+5)03 where 
the two atomic species that alternate with Ti on the B 
site are virtual atoms that we constructed by varying the 
nuclear charge of Ti by ±5. Therefore, as 5 is tuned con- 
tinuously from to 1, we can simulate a set of systems 
evolving from a conventional BaTiC>3 ferroelectric system 
to a heterovalent Ba(Sci/3Ti!/ 3 Vi/3)03 one in which all 



FIG. 2: Energy vs polarization P in Ba(Ti-<5,Ti,Ti+J)0 3 at 
different S. Note the saddle point of E(P) shifts in the direc- 
tion of the secondary (shallower) minimum as 8 increases. 



three alternated species are from neighboring columns in 
the periodic table. 

As a consequence of the compositionally broken inver- 
sion symmetry, the thermodynamic potential associated 
with the FE instability does not have the usual sym- 
metric double-well form. Instead, it takes the form of 
an asymmetric double well, or even of an asymmetric 
single well, depending on the strength of the composi- 
tional perturbation that breaks the symmetry. In nor- 
mal ferroelectrics, it is sufficient to locate one minimum 
of the double-well potential; the other is then obviously 
given by applying the inversion operation. Here, this no 
longer works. A steepest-descent minimization starting 
from the ideal structure typically arrives at the primary 
(deeper) minimum, but the secondary (shallower) mini- 
mum can be rather difficult to find in practice. 

A procedure was described in Ref . |ll| that allows one to 
search for both minima when they coexist. To illustrate 
the procedure, we plot in Fig. [j], for several values of <5, 
the energy as a function of displacement amplitude along 
the straight line in the 15-dimensional parameter space 
connecting the primary and secondary minima. (The di- 
rection along this line is taken to define the "FE direc- 
tion" with £fe = being the midpoint between minima.) 
Unfortunately, however, it is not possible to plot such a 
curve for S > 0.4, since only a single minimum exists in 
this range of S. 

Here, we demonstrate how the current method allows 
for a much more natural treatment of these systems, es- 
pecially at large S. At a fixed value of polarization, we 
calculate F(P) as in Eq. (||). That is, we minimize the 
total energy over the internal coordinates subject to the 
constraint that the spontaneous polarization has a fixed 
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FIG. 3: Calculated polarization at the left minimum (solid 
circle), right minimum (open circle), and saddle point (dia- 
mond) in Ba(Ti-<5,Ti,Ti+5)C>3. Left (right) minimum is the 
principal one for 8 > (8 < 0), as shown in insets. 



value, following the procedure described in Section [I] 



As in our previous work, this is done in a fixed tetrago- 
nal cell. 

The energy as a function of thcpolarization for several 
values of 8 is illustrated in Fig. g. We obtain a similar 
energy evolution as in Fig. [j]. However, there are two im- 
portant qualitative differences. First, the new procedure 
is not limited to the range of S in which both minima ex- 
ist. At larger 8 (e.g., 5 = 0.6), where the secondary mini- 
mum has disappeared due to a strong symmetry-breaking 
perturbation, the new procedure allows the mapping of 
the energy to be carried out just as easily as at smaller 
5. Second, the horizontal axis of the figure now has a 
physical meaning of polarization. For example, a glance 
at Fig. ^ shows an interesting feature, namely that the 
saddle point is also polarized, unlike in a normal ABO3 
compound. 

We investigated this interesting feature further by plot- 
ting in Fig. the polarization at the saddle point, as well 
as at the energy minima, as a function of 8. In the "nor- 
mal" case 8 — 0, the saddle point is unpolarized and 
the two equivalent minima carry equal and opposite po- 
larizations. However, all the stationary points are seen 
to be shifted in the direction of the shallower minimum 
as 8 is turned on. As a consequence of these shifts, the 
polarization of the primary minimum changes sign near 
8 ~ 0.4. At the critical 8 the saddle point and secondary 
minimum meet and annihilate. Returning to small values 
of 8, a closer analysis (not shown) indicates that the po- 
larization at the saddle point increases in proportion to. 
S 3 . This observation agrees well with previous studiesEJ 
showing that certain other measures of the effect of the 
symmetry-breaking perturbation also scale as 8 3 . 

In summary, we have illustrated the convenience and 
power of our new method in the context of recent work 
on a new class of ferroelectric materials with composi- 
tionally broken inversion symmetry. The new approach 
is especially useful for studying the case where the com- 
positional perturbation is so strong that only a single 
local minimum survives. In the former procedure, the 



FIG. 4: Fully relaxed z coordinates for each atom in BaTi03 
(see unit cell at right) as a function of polarization P z in 
the simple-cubic lattice. Curves are cubic-spline fits to cal- 
culated points; top- and bottom-most points correspond to 
translational images of Ba and 03 atoms, respectively, in 
neighboring unit cells. P„ marks the spontaneous polariza- 
tion. The Born effective charge Z* for each atom is marked 
at P z = 1.64C/ni . Shaded area indicates metallic regime. 



definition of the FE direction was based on the location 
of two local minima, and was therefore useless when one 
of the minima had disappeared. On the contrary, in the 
new method the energy surface can be straightforwardly 
mapped out as a function of polarization, regardless of 
whether the secondary minimum exists or not. Moreover, 
expressing the behavior as a function of polarization pro- 
vides a much more informative picture of the system. For 
example, certain interesting and non-trivial behaviors of 
the polarizations at the saddle points and minima can be 
elucidated. 



B. Structural response in BaTiC>3 

In this section, we apply our approach to BaTiC>3, one 
of the most-studied perovskite ferroelectric compounds. 
It undergoes a sequence of structural phase transitions 
with decreasing temperature: from the high-temperature 
cubic to the tetragonal phase at 130°C, then to an or- 
thorhombic phase at 5°C and finally to the ground-state 
rhombohedral structure at — 90°Cx2l The three successive 
phases, distortions of the cubic perovskite structure, are 
characterized by spontaneous polarizations aligned along 
the [001], [011], and [111] directions respectively. 

Here, we focus on the dependence of the internal struc- 
tural parameters of BaTiOa on the polarization P. Our 
calculation is restricted to allow only atomic displace- 
ments along the z direction in a fixed simple cubic lattice, 
with full relaxation of the internal structural parameters 
at fixed polarization P = P z z to yield equilibrium coordi- 
nates R eq (P 2 )- In this way, we can investigate the contri- 
bution of the internal structural parameters alone, decou- 
pled from the strain degrees of freedom, to the structural 
response to an electric field, providing a first step towards 
understanding the nonlinearities of the total structural 
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P 2 (C/rri ) P z (C/m ) 

FIG. 5: The component of the unit displacement vector £ z 
corresponding to each atom in BaTiOa (left) and PbTiO;; 
(right) as a function of polarization P z . In the left panel, 
vertical lines demarcate the regimes of BaTiOs soft-mode- 
like, PbTiOs soft-mode-like, atom-pair, and metallic behav- 
ior; at right, the single vertical line separates soft-mode-like 
and atom-pair regimes. 



response expected with increasing £. 

Experimentally, BaTiOa-is known to have a cubic lat- 
tice constant of 7.547 a.u.Ea Our LDA calculation yields 
an equilibrium lattice constant of 7.45 a.u., 1.3% smaller 
than the experimental value, an error typical of the LDA. 
As is well known, the ferroelectric .instability depends 
sensitively on the crystal volume EJOS We therefore 
choose to work at the experimental cubic lattice constant. 

The spontaneous polarization P s is obtained by full 
relaxation of the internal structural parameters, and is 
found to be 0.21C/m 2 . The relaxed internal coordinates 
for each Ba, Ti, 01 and 03 atom are plotted as a function 
of P z in Fig. [|. For P z > P s , the state can be realized as 
an equilibrium state in an appropriate fixed electric field, 
while states with P z < P s are local maxima of F(P Z ) — 
£P Z for some value of £. For example, the value of £ 
corresponding to P z — 0.48C/m 2 (approximately twice 
P s ) is 16 MV/cm. 

To focus on the dependence of the character of 
the distortion on the amplitude P z , we define a 
"unit displacement vector" (£2 a , £j\ £? 3 ) by nor- 
malizing the sum of the squared displacements to 
one. At P z = P s , the unit displacement vec- 
tor is found to be (0.26,0.73,-0.22,-0.55), closely 
resembling the unstable ferroelectric mode of cu- 
bic BaTi0 3 (0.18,0.74, -0.18,-0.59) computed from a 
linear-response calculation. In Fig. g, we show the P z 
dependence of the components of the unit displacement 
vector. If the polarized state were obtained by freez- 
ing in a single polar mode, these components would be 
constant. The actual behavior is considerably more com- 
plicated. 

Three distinct regimes for the atomic displacement 
pattern can be clearly observed. For P z below w 
0.48C/m 2 , the relative displacements are similar in char- 
acter to those of the soft mode. In this regime, the mag- 
nitudes of the Ba and 01 components increase with P z , 
while the magnitudes of the Ti and 03 components de- 
crease. For P z between roughly 0.48C/m 2 and 1.4C/m 2 , 
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FIG. 6: Computed Born effective charges for each atom in 
BaTiOs as a function of the polarization P z . For the cubic 
structure (-P*=0), we obtain Z^ = 2.72, Z^i = 6.99, Zq ± = 
-5.57, Z* .. = -2.07. 



the consequence of these opposing trends is that the mag- 
nitudes of the Ba and 01 displacements actually exceed 
those of Ti and 03, respectively, changing the character 
of the structural distortion. At P z w 1.4C/m 2 , the trend 
with P z reverses for Ba and Ti, so that as the polarization 
increases further, the Ba and 01 atoms move together in 
a direction opposite to that of Ti and 03. 

The Born effective charges Z* are expected to be sensi- 
tive to the internal structural parameters. Figure]^ shows 
the evolution of the computed Z* for each atom. Near 
the cubic structure, the dependence of Z^ { and Zq 3 on 
P z is nearly quadrati<p-iin agreement with previous cal- 
culations for BaTi03.E?J While Z^ a and Zq 1 are rather 
insensitive to P z , the Born effective charge of Ti decreases 
by over 30% from P z = to 1.4C/m 2 , with a correspond- 
ing increase in that of 03. More specifically, Z^ drops 
to its smallest value +4.7 while the magnitude of Zq 3 is 
close to its smallest value —3.9. This structural sensi- 
tivity can be understood as being related to the anoma- 
lous values in the undistorted cubic perovskite structure, 
which arise from the hybridization of Ti and O orbitals 
in the Ti-03 chains oriented along z. As shown in Fig. |J, 
when P z increases, the Ti are displaced towards one 03 
neighbor and away from the other, disrupting the chain 
and reducing the anomalous displacement-induced cur- 
rent along the chain.a Consequently, the magnitudes of 
the Born effective charges are reduced towards the nom- 
inal valences +4 for Ti and —2 for O. 

As the polarization and associated structural distor- 
tions become larger, the bandstructure evolves corre- 
spondingly. We find significant changes in the band 
structure for P z = 1.8C/m 2 relative to that of the undis- 
torted cubic perovskite structure. The hybridization be- 
tween the Ti 3d and O 2p bands becomes more signif- 
icant, as expected from the decreased Ti-03 distance. 
Some bands, such as the topmost O 2p band, lose the 
characteristie-|flatness that is usually seen in perovskites 
like BaTiOaL 2 ] and become much more dispersive. Most 
importantly, the band gap decreases with increasing P z , 
extrapolating to an insulator-metal transition just above 
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P z = 1.8C/m 2 . The polarization is a meaningful quantity 
only in insulators, and therefore calculation for higher 
values of P z cannot be considered. 

We carried out an analogous calculation for PbTi03 
in the cubic structure using the lattice constant deduced 
from experiment, ao — 7.5 a.u, yielding a spontaneous 
polarization of 0.73C/m 2 . The unit displacement vector 
as a function of P z is shown in the right panel of Fig. ||, 
where the pattern resembles that of BaTiC>3 at interme- 
diate values of P z . Thus, the pattern that is field-induced 
in BaTiC>3 is characteristic of that of PbTiC>3 at zero elec- 
tric field. The high-P z regime sets in at around 1.5C/m 2 , 
with the Pb and 01 atoms moving together as a pair and 
Ti and 03 moving together as a second pair. 

Some general observations can be made about the ef- 
fects of an electric field on the internal structural param- 
eters. At small fields, the cations and anions move inde- 
pendently, following the electrostatic force correspond- 
ing to the sign of the formal valence. However, once the 
field-induced distortions are large enough so that short 
range interatomic repulsion prevents further compression 
of Ba/Pb-01 and Ti-03 bonds, the further distortions 
acquire a new character in which these atoms move as 
pairs. From the computed values of the dynamic effec- 
tive charges at P z = 1.64C/m 2 (see Fig. gj), the Ba-01 
pair and Ti-03 pair carry net charges very close to —2 
and +2 respectively. 

In summary, we have shown that the "simple" per- 
ovskite compound BaTi03 exhibits significant nonlinear- 
ity in structure with increasing polarization, correspond- 
ing to large electric fields, while the atomic displacement 
pattern of PbTi03 is relatively slowly varying. It is in- 
teresting to note that at large enough polarization, the 
atomic displacement pattern of BaTi03 in fact resem- 
bles that of PbTi03. This may help to shed some light 
on the factors responsible for the differences in properties 
between alkaline earth and Pb-based perovskites. 



C. Non-Linear Dielectric and Piezoelectric 
Response 

Tunability of the dielectric and piezoelectric coeffi- 
cients by an applied electric field, a property of great 
technological importance, is expected to be especially 
large in ferroelectrics due to the dependence of these coef- 
ficients on electric-field-induced structural changes such 
as those reported for BaTi03 in the previous section. 
This behavior can be quantified by the values of the non- 
linear dielectric and piezoelectric coefficients. In this sec- 
tion, we formulate the calculation of these nonlinear co- 
efficients in our polarization-based framework, and give 
results for tetragonal PbTi03. 

The first step in this analysis is the computation of 
F(P) and fj(P) from the minimization of F(R, rj, P) 
at fixed P. This is followed by the minimization of 
F(P) — P £ at fixed £, directly yielding P(£) and 
ij(£) = rj(P(£)). From the first derivative of P(£ ), we 



obtain the field-dependent static dielectric susceptibil- 
ity tensor Xa@{£ )) with the nonlinear coefficients defined 
through a small-£ expansion 



1 dP a (£) 



= X$ + X$ 7 £r+0(n- ( 31 ) 



,(2) 



i)£ 



The relative dielectric tensor is given by e a p — S a /3 + 
Xa/3- Correspondingly, from the first derivative of rj(£), 
we obtain the field-dependent piezoelectric tensor dy,p{£), 
with the nonlinear coefficients defined through a small-£ 
expansion 



d£ ^ 



d {2) £ 



0{£ 2 



(32) 



In fact, in our present implementation we perform the 
minimization of F(R, rj, P) at fixed P in two separate 
steps. First, we obtain a reduced free-energy density 
F(rj, P) by minimizing with respect to R at fixed r\ and 
P. Further minimization with respect to 77 to obtain 
F(P) and i](P) allows the computation of the zero-stress 
responses as in the previous paragraph. In addition, this 
approach allows the computation of the clamped-strain 
dielectric response, measured at frequencies above the 
resonant frequency of the sample, through minimization 
of F(rj, P) — P • £ at fixed £ and rj — rj(£), directly yield- 
ing P{rj{£),£) and x(v(£)i£)- This two-s tep pr ocedure 
is also required, as mentioned in Section III A , by the 
limitations imposed on the present implementation by 
the use of ABINIT 3.1. For present practical purposes, 
F(r), P) is obtained in a parameterized form by fitting a 
Landau-Devonshire expansion to values of F obtained by 
calculations for an appropriate set of r\ and P. 

We have applied this procedure to compute the non- 
linear dielectric and piezoelectric response of tetrago- 
nal PbTi03 to fields along z using the i = 1 expres- 
sions (Eq. ^). At this level of approximation, only lat- 
tice contributions to Xap{£) and dy.p(£) are included, 
and their electric-field dependence arise only through in- 
duced structural changes. However, this is expected to 
be a good approximation for PbTi03, where the lattice 
contribution to the dielectric and piezoelectric responses 
dominates even at T = 0. 

With a field along z, symmetry constrains the struc- 
tural response to consist of a tetragonal strain, specified 
by two independent parameters r/i = r\ xx = r\ yy and 773 = 
rj zz , and a set of atomic displacements along the I di- 
rection described by three independent parameters. Cor- 
respondingly, we obtain an expression for F(R, 77, P) = 
F(r/i,ri3, P z ) by computing minR E(R, rj, £ = 0) for a set 
of tetragonal cells with the constraint P(R, 771, 773) = P z z. 
The results are used to fit the parameters in a Landau- 
Devonshire expression expanded around the minimum- 
energy cubic structure (ao=7.33 a.u.), 



F(m,V3,Pz 



= E + -d^nl 



A-20oP z + ^400 P z "I - ^600 

ZBiyyrjiP* + B lzz rj 3 P* , 



) +(712(277! 773 
AfinnPf 



(33) 
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TABLE I: The values of the least-squares fitted parameters in 
Eq. (^3|) at £ = in PbTi03. The units are the appropriate 
combinations of Ha and(C/m 2 ) 2 . 

Parameters Values Parameters Values 



E 


-165.953 


A400 


0.005 


Cxi 


4.374 


^600 


0.004 


Cl2 


1.326 


B\ zz 


-0.199 


A200 


-0.003 


B\yy 


-0.049 



TABLE II: Comparison of the structural parameters com- 
puted by minimizing Eq. (^3|) with those computed 
from direct— iHDA calculation and those obtained from 
exp eriment OE3 

AE (mHa) P/=° (C/m 2 ) a (Bohr) c (Bohr) 



model 
LDA 

exp 



0.86 
0.90 



0.67 
0.65 
0.75(295K) 



7.324 7.487 
7.310 7.484 
7.373 7.852 



where P z is the polarization per unit volume and the 
truncations to sixth order in P z and to lowest order in 
the elastic and polarization-strain coupling are found to 
be sufficient within a standard least-squares fit. The re- 
sulting coefficients are shown in Table ||; statistical anal- 
ysis shows that the strain coupling parameters B\ yy and 
Bizz are the most sensitive to changes in the input con- 
figuration energies. 

We now use this expansion to compute the field depen- 
dence of the strain and polarization under zero stress by 
minimizing F(rji, 773, P z ) — £ -P z with respect to 771, 773 and 
P z to get 771(f), 773(f) and P z {£)- By first considering £ 
= 0, we can confirm the validity of the parameterization 
by comparing the tetragonal structure obtained by mini- 
mizing the expression for F (771, 773, P z ) given by Eq. J33] ) 
with properties of the fully relaxed tetragona l g round- 
state structure in zero electric field. In Table [n] we list 
the energy difference AE between the tetragonal ground 
state and the cubic structure, the spontaneous polariza- 
tion P z (£ = 0), and the lattice parameters, finding good 
agreement in all respects. 

Next, we consider nonzero £. Minimizing first with 
respect to 77 gives a free energy 



F{P Z ) = A 2m P 2 z + A i00 P* + A eao P? 



(34) 



where 



A400 = A. 



1c\iB\ zz B\ vv — c\\B\ — |( c n + C \2)B\ 



400 



(cn + 2ci 2 )(cii 



C12) 



(35) 



and A400 is found to be 4.5x 10~ 4 HaC~ 4 m~ 8 . 

Since A 2 oo < 0, F(P Z ) has a double well structure, so 
that F{P Z ) — P z £ has two local minima for small enough 
values of £. The evolution of the two local minima with 
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FIG. 7: Calculated polarization-vs.-electric-field hysteresis 
loop (upper panel) and static susceptibility (lower 
panel) of PbTiOa under stress-free condition (solid circle) and 
clamped-strain condition (open square). Dashed line corre- 
sponds to the non-accessible state (saddle point in the ther- 
modynamic potential). 



TABLE III: Comparison between theory and experimenltll 
(at room temperature) for the first- and second-order dielec- 
tric constants of PbTi03 (RT). The superscripts a and r\ in- 
dicate whether the measurement is under constant-stress or 
constant-strain condition. 

X33 X V 33 xft (nm/V) xiT(nm/V) 



model 67 37 
experiment 79 33 



315 



82 



£ can be summarized in the calculated hysteresis loop 
shown in the upper panel of Fig. [?]. We find an intrinsic 
coercive field £ c of 1.5MV/cm. From Eq. (j3l|), we can 
proceed to calculate the static susceptibility X33(£) an d 
the result is plotted in the lower panel of Fig. ^ Fitting 
this to Eq. (|3l|), we find that the zero-field stress- free sus- 
ceptibility X33 is X33 = 67, the superscript a indicating 
stress-free conditions. 

For the clamped-strain response at zero field, we fix 
77 at rj(£ — 0). A different double well structure is ob- 
tained for F™ (P z ), resulting in a different hysteresis loop 
shown in the same figure. We find an intrinsic coercive 
field £ c of 3MV/cm. From fitting to Eq. (|l|), we obtain 
t^'j = 37, the superscript 77 indicating the clamped-strain 
condition. In Table III , these values are compared with 
the reported experimental dielectric constants i-ajb both 
the constant stress and clamped-strain conditionEH which 
were measured below and above the sample resonant fre- 
quencies respectively. The value for X33 can also be, com- 
pared with a previous first-principles calculation.E3 

In both the free-stress and fixed-strain case, the hys- 
teresis profile of the static susceptibility shows that %33 
increases with field amplitude for the local minimum at 
£ < £ c and decreases with increasing field for the global 
minimum, which is the only branch in the region above 
£ c . For each branch, we find a non-linear susceptibil- 
ity X33 of magnitude 315nm/V in the stress-free case. 
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FIG. 8: The calculated equilibrium strains 771 and 773 and the 
piezoelectric tensor di3 and CZ33 as a function of the electric 
field in PbTiQ 3 . 



However, when the strain is clamped, the coercive field 
becomes larger than in the stress-free case, and the non- 
linear susceptibility is more than two times smaller. In 
the present framework, this is not surprising since the 
change in the dielectric response is the result of a field- 
induced change in structure, and this change is reduced 
by clamping the strain. In nonzero field, the suscep- 
tibility can be either of these two values depending on 
whether the the system is in a single domain correspond- 
ing to the global minimum or to the local minimum, or an 
intermediate value if both types of domains are present. 



D. Field-Induced Structural Phase Transitions 



In a single crystal, the relative stability of distorted- 
structure phases with polarizations in different directions 
is expected to change as an electric field is applied. In 
particular, a phase transition might be induced by ap- 
plying a sufficiently large field in a different direction 
from the polarization of the ground state. This change in 
phase, the result of an electric-field induced rotation of 
the polarization, may be accompanied by a large change 
in strain, manifested as a large piezoelectric response. 

This "polarization rotation" mechanism was proposed 
in Ref. |l3| to explain the experimentally observed colos- 
sal piezoelectric response to electric fields along [001] 
of single-crystal rhombohcdral perovskite alloys such 
as [Pb(Zn 1/3 Nb 2/ 3)03](i_ K )-[PbTi0 3 ] ;c (PZN-PT) with 
compositions near the rhombohedral-tetragonal mor- 
photropic phase boundary (R-T MPBj) , and has been the. 
subject of continuing experimentaEJ and theoreticalEll 
investigation. Particular attention has focused on the 
nature of the path followed by the polarization vector 
with increasing field strength. An effective Hatailtonian 
study of PbZr^Tii^Os near the R-T MPBEl showed 
that with increasing electric field along [111], the po- 
larization vector of tetragonal PZT rotates continuously 
from the tetragonal [001] direction to the rhombohedral 
[111] direction through a monoclinic "Ma" phasec3 with 
P along [mtl]. In contrast, for the case of an [001] electric 
field applied to rhombohedral PZT, the polarization vec- 
tor does not simply follow the return path, but instead 
follows a-discontinuous path of a kind first discussed by 
Nohcda.c3 It first rotates continuously into the Ma phase 
for small field strengths, and-then jumps discontinuously 
to a monoclinic "Mc" phasala with P along [uOl] before 
reaching the tetragonal structure. The calculations show 
that a large piezoelectric response is expected for this 
latter type of path. 



Next, we consider the piezoelectric response (Eq. 32). 
In Fig. ||, we plot the equilibrium values of the strains 771 
and T73 as a function of the electric field along the z di- 
rection. The slopes of these curves give rise to the piezo- 
electric coefficients di3 and g?3 3 which are plotted in the 
lower panels of the same plot. We find e? i3 = — 0.6pC/N 
and ^33 = 40pC/N, considerably less than the room tem- 
perature values measured experimentallyLj (— 25pC/N 
and 117pC/N, respectively) and computed from first 
principlesO We attribute this primarily to the choice of 
pseudopotentials, which give a low value for the ground 
state tetragonal ratio c/a and in particular, a value of a 
almost unchanged from the cubic ciq. However, our calcu- 
lation does serve to demonstrate the applicability of our 
method to the calculation of these quantities. In partic- 
ular, there is to our knowledge no previous calculation of 
the nonlinear piezoelectric response. 



In this section, we apply the full three-dimensional for- 
malism described in Section II to study, in a Pb-based 
perovskite system, the rotation of polarization by an ap- 
plied electric field in the two cases most relevant to en- 
hanced piezoelectric response near the R-T MPB: (i) ap- 
plication of an electric field along [111] to a tetragonal 
system, and (ii) application of an electric field along [001] 
to a rhombohedral system. For (i), we consider tetrag- 
onal PbTiC>3. For (ii), we introduce a simple modifica- 
tion of the structural energetics of PbTiC>3 to stabilize 
a rhombohedral ground-state structure. This follows the 
spirit of a view of PZN-PT and PMN-PT as large-strain 
PbTiC>3-based systems that have been chemically "engi- 
neered" to make them marginally stable in the rhombo- 
hedral phase. CJ We do something very similar, but using 
a theoretical manipulation that avoids the unnecessary 
complexities of the real alloy systems. 
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Free-energy functional 



IV C 



uate 



Extending the procedure described in Section 
to the full three-dimensional case, we first eva 
the reduced free-energy function P(7/,P) by minimizing 
P(R, ?/, P) with respect to R for a set of selected tetrago- 



nal, rhombohedral and orthorhombic structures. Strains 
arc defined relative to the cubic structure with the ex- 
perimental lattice constant (ao = 7.5 a.u.). In the range 
of rj and P of interest, we used a procedure similar to 
Sec. IV C to fit F(rj 1 P) in a Landau-Devonshire form: 



F(r),P) = E + Cx(r,x + m + m) + ^C n (r] 2 + + vt) + C 12 {r, 2m + VsVi + V1V2) + ^f^l + vt + vl) 

+ A 2m {P 2 x + Pi + P x 2 ) + A 400 (Px +Py+ P*) + M^(P 2 y Pl + if Pj + PlPl) 

+ A 60Q (P! + Pi + Pf) + A 420 [P 2 (P£ + Pt) + Py(P* + P*) + P 2 *{P$ + Py)} + M22P 2 x P 2 y P 2 z 

+ B lxx (viP 2 + mP 2 y + nzPl) + B lyy [ m (P* + p 2 z ) + m {p 2 z + Pi) + m {p 2 x + p y 2 )] 

+ B 4yz ( m P y P z +r, 5 P z P x +7 l6 P x P y ). (36) 



We list the least-squares fitted coefficients in the column 



denoted by Ml in Table [IV 
Using Eq. (| 



3|), we now consider the energetics of states 
with different orientations of the polarization in zero 
field. Specifically, we consider E(9, 0, £ = 0), obtained 
by fixing the direction of P along the direction specified 
by spherical angles 9 and relative to the polar axis z, 
and minimizing P(ry,P) with respect to the strain and 
the magnitude of P. As shown in Fig. ^, the tetrago- 
nal phase (T) with polarization along [001] is the global 
minimum, with a saddle point at the orthorhombic phase 
(O) with P || [110] and a maximum at the rhombohedral 
phase (R), with P || [111]. As shown in Table [v|, the 
structural parameters and the spontaneous polarizations 
agree well with the LDA results, especially for the O and 
R phases. 

From this table, it can also be seen that the energy dif- 
ferences between the T, O and R phases are quite small. 
For this reason, the parameters obtained by a global 
least-squares minimization do not accurately reproduce 
the LDA values. In particular, the energy difference be- 
tween the T and R phases is seen to be much larger than 



TABLE IV: Least-squares fitted values of the parameters in 
Eq. (pii]). Ml, all parameters freely varied; M2, with the con- 
straint A222 = 0.062Ha C" 6 m" 12 (boldface). Units are the 
appropriate combinations of Ha and (C/m 2 ) 2 . 





Ml 


M2 




Ml 


M2 


E 




165.947 


-165.947 


Ci 


0.168 


0.168 


^200 




-0.01 


-0.009 


Cn 


3.829 


3.973 


^400 




0.008 


0.005 


C12 


1.462 


1.484 


A22O 




0.015 


-0.0007 


C44 


1.174 


1.218 


^600 




0.003 


0.004 


Blxx 


-0.235 


-0.234 


A42O 




0.010 


0.019 


Blyy 


-0.048 


-0.0525 


A222 




0.009 


0.062 


B^yz 


-0.069 


-0.068 



the LDA result. However, these features of the energy 
surface are crucial to the physics of the structural phase 
transitions. Therefore, we adjusted the fitting procedure 
to reproduce these relative energies accurately while us- 
ing the least-squares procedure for the best overall fit to 
the remaining data, as follows. Rather than introduce 
additional parameters by including higher-order terms, 
we "tune" one parameter while determining the other 
13 parameters by standard least-squares minimization, 
and choose the value for the single tuned parameter that 
yields the most accurate values for both the O-T and O- 
R energy differences. A 222 proves to be the best choice 
for the tuning parameter, and with ^222= 0.062 and the 
other parameters as given in the column denoted by M2 
of Table IV, both the O-T and O-R energy differences 
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FIG. 9: Contour plot oiE{6,4>,£ = 0) for PbTiOs. Spherical 
angles 6 (polar) and <f> (azimuthal) indicate the direction of 
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TABLE V: Comparison between the structural properties of 
the tetragonal (T), orthorhombic (O), and rhombohedral (R) 
phases of PbTiOa. LDA denotes direct LDA structural relax- 
ations; Ml and M2 are as in Table [F^. Units of polarization 
P, rhombohedral angle a, cell volume V, and phase energies 
E axe C/m 2 , de grees, Bohr 3 , and mHa, respectively. 





LDA 


Ml 


M2 


Vt 


399.9 


402.3 


401.9 


c/(l 


i no/i 

1.UZ4 




1 HQ 

l.Uo 


Pt 


0.65 


0.75 


0.71 


Vr 


398.4 


397.3 


398.7 


an 


89.7 


89.6 


89.8 


Ph. 


0.33 


0.32 


0.34 


Vo 


398.8 


398.1 


399.7 


ao 


89.5 


89.4 


89.6 


Po 


0.41 


0.42 


0.44 


En-Eo 


0.060 


0.097 


0.064 


Eq — Et 


0.159 


0.639 


0.154 



as well as the structural parameters and spontaneous po- 
larizations of all three phases are in excellent agreement 
with the LDA results, as shown in the last column of Ta- 
ble [y| Therefore, this set of parameters was used in the 
following calculations. 



2. Engineering a Rhombohedral Structure for PbTiOa 

In previous first-principles investigations of PbTiOs, it 
has been observed that the strain coupling is responsible 
for stabilizing the tetragonal ground state structure.cfl 
In the simple cubic lattice, the lowest-energy structure 
has polarization along [111], corresponding to a rhombo- 
hedral symmetry, while the energy of the optimal state 
with polarization along [001] is higher. However, when 
the lattice is allowed to relax, the energy gain from strain 
coupling in the tetragonal structure is much larger than 
the gain in the rhombohedral structure, leading to the 
observed reversal of stability. In each case, the energy 
gain from strain coupling increases as the relevant elastic 
constant decreases. So, if it were possible to decrease the 
shear modulus C44, there would be a critical value below 
which the rhombohedral state would be most stable. 

Within Eq. (|36|), the modification of C44 can be imple- 
mented by the inclusion of a tunable shear elastic term 

F( V ,P) = F( V ,P) + ±AC M (vl +Vl +Vl) (37) 

where AC44 — corresponds to PbTiC>3 with its natu- 
ral shear elastic modulus. Using Eq. (|37|), we compute 
the zero-field energy for the optimal tetragonal, rhombo- 
hedral and orthorhombic phases as a function of AC44. 
This yields the phase sequence shown in Fig. [l^ with the 
T and R phases separated by a sliver of an orthorhombic 
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FIG. 10: The energies of the rhombohedral (square) and 
tetragonal (circle) phases relative to the orthorhombic phase 
(chosen as the zero of energy) as a function of the tunable 
shear modulus AC44 of PbTiOs, calculated using Eq. (fjTj). 
The ranges of AC44 in which the tetragonal, orthorhombic 
and rhombohedral phases are most stable are separated by 
vertical dashed lines and indicated by T, O and R respec- 
tively. 
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FIG. 11: Same as Fig. hoj but calculated using a direct LDA 
approach. As in Fig. |10^ an orthorhombic window appears, 
though the phase boundaries are slightly shifted. 

phase. This phase sequence is very reminiscent of that of 
the Pb(Zn 1/3 Nb2/3)0 3 -PbTi03 systerr£3 with the tun- 
able parameter being the proportion of PbTiC>3 . The sta- 
bility of the orthorhombic phase reflects the importance 
of the sixth-order terms in Eq. (|36|), as in a fourth-order 
model only tetragonal and rhombohedral structures are 
possible minima. 

To check that the observed phase sequence is not an 
artifact of our fit, we have computed the structural pa- 
rameters and energies of the tetragonal, orthorhombic 
and rhombohedral phases as a function of AC44 through 
direct LDA calculations. For consistency with Eq. (|3^), 
we implement the adjustment of the shear modulus as an 
additional applied stress 

Act, = -AC 44 r?i (38) 

where ui (with i = 4, 5, 6) are the shear stress compo- 
nents in Voigt notation. The results, given in Fig. |ll| , 
show the same T-O-R phase sequence as Fig. [n| While 
the T-0 and O-R phase boundaries are slightly shifted, 
the width of the orthorhombic window is comparable to 
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FIG. 12: Contour maps of E(8, <f>, £) on the upper hemisphere 
< 9 < ^tt for an electric field of magnitude £ applied along 
the pseudocubic [111] direction to tetragonal PbTi03. The 
contour is equally spaced in log(E — -E m m + 8), where -E m i n is 
the global minimum and 8 is a small offset. The central axis 
points along the [111] direction, (a)-(d) correspond to electric 
fields of 0, 0.86, 1.73, and 3.46xl0 3 kV/cm, respectively. 

that in Fig. [H]. Thus, in the following, using Eq. ( pTIj ) 
with a particular value of AC44 , we expect results which 
would reflect a direct LDA calculation, though perhaps 
with a slightly different AC44. 

3. Electric-field-induced phase transitions 

In a single crystal, the relative stability of phases with 
polarizations in different directions is expected to change 
as an electric field is applied. In particular, a phase tran- 
sition might be induced by applying a sufficiently large 
field in a different direction from the polarization of the 
ground state. Here, we consider two such cases: tetrago- 
nal PbTiC>3 in an electric field along [111], and rhombo- 
hedral "PbTi03," stabilized by a nonzero value of AC44, 
in an electric field along [001]. 

First, we consider tetragonal PbTiC>3 (AC44 = 0) in 
an electric field along [111], which tends to favor a rhom- 
bohcdral direction for the polarization. To investigate 
the evolution of various phases with fm, where fm is 
the magnitude of the electric field, we perform the mini- 
mization in two steps. First, we transform the Euclidean 
coordinates (P x ,P y ,P z ) into spherical coordinates (P,8,<p) 
and compute 

E(6, <t>, £ m ) = mm[F(v, p ) - ( p x + Py + Pz)/Vs] 

P,rj 

(39) 
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Electric field (10 3 Kv/cm) 

FIG. 13: The Cartesian component P x (circle), P y (triangle), 
P z (square) of the polarization as a function of the magni- 
tude of the electric field applied along the [111] pseudocubic 
direction in PbTiOa. Inset shows the polarization path. 



Then, we locate the minima on the sphere of polarization 
directions parametrized by 9 and <p. 

The evolution of the phase stability can be readily 
displayed by the contour plots of E{9, <f>, <?m) shown in 
Fig. [lj. At zero electric field, the tetragonal structure 
appears as a three-fold degenerate energy minimum in 
the hemisphere shown. As £111 increases, the minima 
migrate from the tetragonal positions along the lines cor- 
responding to the monoclinic Ma phase (three-fold de- 
generate) and eventually reach the rhombohedral point 
at the center of the hemisphere. 

Figure ^ shows how the polarization components of 
tetragonal PbTiC>3 evolve with the amplitude of £\u. 
At fin = 0, the only non-zero component is P z . As £ 
increases, P x = P y grow while P z slowly decreases. The 
structure thus enters the Ma monoclinic phase. When 
£ reaches 1.4 x 10 3 kV/cm, the three components merge 
and the system enters the rhombohedral phase where the 
polarization vector points along the pseudocubic [111] di- 
rection. While rotating, the polarization vector remains 
in the (110) plane, as shown in the inset of Fig. |l3| . 

Next, we consider rhombohedral "PbTiOs" with 
AC44 = — l.IHa (see Fig. |l(]) in an electric field along 
the [001] direction, which tends to favor a tetragonal di- 
rection for the polarization. The analogue of Eq.( ||^) 
is 

E(0, 0, £001) = mia[F{r), P) - £ Q01 P z ] (40) 

P,TJ 

The energy contour plot in this case is shown in Fig. |lj. 
In zero electric field, the system is in a rhombohedral 
phase with an eight-fold degenerate minimum. For small 
nonzero foot, the energy minima correspond to a Ma 
phase as shown in Fig. |l4j(b-c) where there are four de- 
generate minima lying in the (110) plane. At a critical 
value of £001 ) the energy minima jump to four- fold points 
in the (100) plane, as can be seen in Fig. |l4|(d). The four 
minima then move smoothly towards the [001] axis, fi- 
nally merging to yield the tetragonal phase. 
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FIG. 14: Contour map of the free energy (upper hemisphere) 
when an electric field is applied to rhombohedral "PbTiOa" 
(AC44 = — l.IHa) along the pseudocubic [001] direction. The 
central axis corresponds to one of the tetragonal directions. 
(a)-(f) correspond to electric field magnitudes of 0, 1.4, 2.8, 
7, 14, and 19 x 10 3 kV/cm, respectively. 



Figure [15| shows how the polarization components of 
rhombohedral "PbTiOs" evolve with the amplitude of 
£001 ■ Under zero applied electric field, the polarization 
vector starts along the pseudocubic [111] direction (P x = 
P y = P z = 0.56C/m 2 ). As £ oi increases, the structure 
enters an Ma phase in which P x and P y remain equal, but 
become less than P z . P x and P y keep dropping until P y 
shows a sudden jump to zero at around 4.5 x 10 3 kV/cm. 
At the same time, both P x and P z exhibit an upward 
jump in their values. The new phase corresponds to a 
different monoclinic phase denoted by Mq- The structure 
remains in the Mq phase until P x also drops to zero 
at around 19 x 10 3 kV/cm, yielding a tetragonal phase. 
As the field increases further, P z continues to increase 
smoothly. 

In this section, we have seen that a small modification 
of the structural energetics of PbTiC>3 can yield a com- 
plex polarization path quite similar to that proposed by 
Nohedao and observed in simulations of PZT.cil Addi- 
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FIG. 15: Same as Fig. [l^, for an electric field applied along 
the [001] pseudocubic direction in rhombohedral "PbTiOV 
obtained with AC44 = —l.IHa. 

tional calculations, for example of the lattice parameters 
as a function of electric field, may assist in achieving a 
direct experimental observation of this behavior. In ad- 
dition, further exploration within this framework may 
suggest ways to produce and control complicated polar- 
ization paths in real systems. 



V. SUMMARY 

In this paper, we have introduced a formalism for com- 
puting the structural response of an insulating system to 
a static homogeneous macroscopic electric field. We have 
shown that, in the presence of an electric field, the ther- 
modynamic potential E(R, 77, £ ) can be minimized by in- 
troducing a related thermodynamic potential F(R, 77, P) 
in which the polarization P is treated as a fundamen- 
tal variable. Corresponding to each polarization P, the 
equilibrium values for the internal coordinates R and 77 
as well as the minimum of this energy functional can be 
computed. Consequently, one arrives at an energy func- 
tional that only depends on P and where the effect of 
a homogeneous electric field can be treated exactly by 
adding a linear term — £ ■ P to this functional. 

In practice, when E(R, 77, £) is expanded to first order 
in £, the minimization is reduced to one over the internal 
coordinates constrained by a fixed polarization computed 
at zero electric field. We have implemented a minimiza- 
tion scheme in the framework of a modified DFPT, using 
a consistent discretization formula that was developed 
for the response to an electric field. Consequently, the 
computed response is compatible with the Berry-phase 
polarization, which is a central quantity in the formal- 
ism. 

It is important to note that the present 7 = 1 theory 
is most useful for systems in which the response to an 
electric field is dominated by the changes in atomic co- 
ordinates and strains rather than by electronic polariza- 
tion. Ferroelectric and nearly-ferroelectric materials are 
among the best examples. We therefore look forward to 
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future applications of our new approach for a variety of 
purposes, for example, ferroelectric alloys and ferroelec- 
tric superlattices. Applying the method to the so called 
"high-K materials" to study their dielectric properties in 
the presence of an applied electric field also appears to 
be a promising direction. 

Though the higher-order (say i = 2) theory requires 
higher (> 3) order energy derivatives, this does not pre- 
clude its application. As mentioned in Appendix [A|, it is 
possible to approximate certain response quantities that 
are related to the third derivatives by constant values 
from a single structure, if they show only small varia- 
tions within the range of the polarization studies. Sys- 
tems that may satisfy such a condition will be the subject 
of further investigation. 



Acknowledgments 

This work was supported by ONR Grants N00 14-97- 
1-0048 and N00014-00-1-0261. The work of K.M.R. was 
performed in part at the Aspen Center for Physics. We 
would like to thank X. Gonze for his interest in the work 



and valuable discussions on the Abinit code. We ac- 
knowledge Ph. Ghosez and M. Veithen for their help on 
the FHI pseudopotentials. 



APPENDIX A: SECOND-ORDER EXPANSION 
FORMALISM 

This appendix presents the formalism in Section II 
for truncation of the sum in Eq. (^) at i = 2, that is, 
at second order in the electric field £. At this order, 
the thermodynamic potential E(R,r/,£) is replaced by 
E%(R, rj,£), which is the sum of the first three terms in 
Eq. (|). 

We recall the definition of the dielectric susceptibility 
tensor 



1 d 2 E(R,i 1 ,£) _ 1 dP a {R,T],£) 
e d£ a d£rj 



d£ a 



(Al) 



Therefore, we can write 



J 



E 2 (R, 77, £) = E(R, r),£)-Yl Pq ( R ' ^ - f £ £a£ /3X«/3( R ' '?< °) 

a a/3 



(A2) 



and 



P 2m (R, 77, £) = P a (R, 77, 0) + e ^ E 0X a p (R, V, 0) • 





(A3) 



The computation of F(P) (Eq. g) for a given P proceeds by the minimization of E(R, 77, A) + A • P following the 
procedure in Section IIIA. This involves computing the derivatives 



8E 2 (R, 77, A) _ &E(R, 77, 0) ^ dP a (R, 77, 0) 



dR 



dR 



E 



OR 
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^-7^ A Q A ( 



<9Xc«/3(R,z7,0) 



a/3 



dR 



(A4) 



dE 2 (R, v ,\) 5£(R,77,0) 



dP a (R, v ,0) 



A, 



e ° x^ \ \ ^XqMR-^o) 

1 a/3 



(A5) 



9£; 2 (R,77,A) 



-P a (R, 77, 0) - e 2^ Xa/j(R. V, 0)A/3 + P a 





(A6) 



r 



These are related to the corresponding derivatives in 
the 7 = 1 case (Eq. [n]) by the add ition of terms one 
order higher in A. From Eq. (A6), we see that at 
this order P(R, 77, A) includes an electronic contribution 
e o Xa/3(R, 0)A/3. The effective forces and stresses 
(Eqs. lAl and |A5|) involve the derivatives of \ with re- 



spect to R and 77. While these are in principle obtainable 
from the 2n+l theorem, they are not routinely calculated 
in current DFPT codes. For cases where the lattice con- 
tribution to P dominates, it is reasonable however to ap- 
proximate the R and 77 dependence of x by evaluating it 
at the zero-field equilibrium structure. A more accurate 
but still practical approximation would include the first 
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order changes with respect to <5R and 5rj, with the deriva- 
tives computed through a finite difference approach. 



APPENDIX B: MULTIBAND DISCRETIZATION 
FORMULA 



In Sec. IIIB we presented a finite- difference formula, 
Eq. (|30|), representing the derivative id\uk)/dk in the 
single-band ID case. In this Appendix we generalize the 
derivation in order to obtain a corresponding formula for 
the multiband 3D case. 

The general expression for the electronic polarization 
in 3D is easily reduced to a sum of ID Berry phases over 
strings of k points.E3 We can write 



1 



VN k 



(Bl) 



where V is the cell volume, a labels the three primitive 
real-space lattice vectors R a conjugate to the primitive 
reciprocal-space vectors G Q , and kj_ runs over a 2D mesh 
of iVfc positions in the reciprocal-space directions perpen- 
dicular to a. The contribution from the string S(fex) of 
k-points running parallel to G a at a given k± is 



P a (kj 



2tt 



^2 ImlndetM( k < k+b ), (B2) 



k£S(k ± ) 



where / = 2 for spin, 



^mn k+b) = («mk|«n,k+b> 



(B3) 



is the overlap matrix formed of inner products between 
Bloch orbitals on neighboring k-points on the string, b is 
the separation between neighboring points on the string, 
and m an d n run over the occupied valence bands. Equa- 
tion (|B2| ) is ess ential ly the multi-band generalization of 
Eq. (||) of Sec. [hII. 



For the remainder of this Appendix, we drop the 3D 
notation and start from the ID version 



P = -— Vim In det M ( 

9tt Z— I 



k,k+b) 



27T 



(B4) 



of Eq. (jB2j), and correspondingly for Eq. (B3). Our task 
is to compute the variation SP arising fro m th e first-order 
variations of the wavefunctions in Eq. (B4). Focusing 
on a single wavevector k and its neighbor k' = k + b 
and letting M = M^ k ' k \ our central task is clearly to 
compute the first-order variation of the phase 



Im In det M . 



Using 



detM = £(-r/n 



(B5) 



(B6) 



where p runs over all possible permutations among the 
occupied bands, the change in this phase from a first- 
order change in the wavefunctions at k is 



= Im- 



8 dct M 



(B7) 



det M 
where 

SdetM = V^(-l) p V^(<5u„ fc |up ( „) fc ,) {u m k\up{m)k>) 

p n m^n 

= V^(-l) P V^(<5u„ fc |w p(n)fe /) ^0,mp(m(B8) 
p n m^n 

Here Mq is the matrix M evaluated before variation of 
the wavefunctions. 



Unfortunately, Eq. ( |B8| ) does not lend itself to simple 
evaluation. However, we can reduce Eq. (B8) to a trivial 
form as follows. Consider a linear transformation 



\u nk > 



(B9) 



among the occupied states at k', where A is a non- 
singular (but not necessarily unitary) matrix. Letting 

Mmn = (u m k\u n k'), it follows that M = MA and thus 
det(M) = det(M) det (A). Since A is a constant matrix, 



S In dct M = 6 In det M . 



(BIO) 



We thus have the freedom to evaluate Eqs. (B7) and ( ]B8| ) 
with the substitutions M — » M, M — > Mq and u m k> — > 
Umk', where M = MA and Mq = MqA, for arbitrary A. 
The obvious choice is A = M^ 1 . We then find that the 



only permutation that surviv es in Eq. (B8) is the identity 
and the denominator of Eq. (fB7|) becomes unity, so that 



where 



\U n k> 



{Mq ^)mn\Umk' 



Eq. (Bll) can also be written neatly as 
5<j) = ImTr(5M ■ M^ 1 ) . 



(Bll) 



(B12) 



(B13) 



Carrying out similar manipulations for the connection 
between k and k — b, we can define 



\Vnk) 



I 

2f> 



{\u n ,k+b) - \UnM-b)) 



(B14) 



which becomes the finite-difference representation of 
id\u n k) / dk in the multiband case, analogous to Eq. (|3C|) . 
It is easy to check the orthogonality of the v n k to the 
occupied subspace, 



{Unk\V mk ) 



Snm) = , 



(B15) 
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thus removing the need for explicit application of a 
conduction-band projector onto the \v nk ) when comput- 
ing the right-hand side of Eq. (p5|). Since (u m k \u n k' ) = 
5mm we can think of \u n k') defined in Eq. (B12) as 



a phase-aligned and amplitude-corrected "partner" to 
\u n k) formed from the occupied subspace at k' , and \v n k) 
is proportional to the difference between the "partners" 
at k + b and k — b. 



Finally, the variation of Eq. (B4) becomes 



SP = — V* Re (Su nk \v nk 

7T * — ' 

k 



(B16) 



in analogy with Eq. (J2^) . 

Our implementation of this scheme into ABINIT is 



based on Eqs. (B11-B16) above. 
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